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Abstract 

We investigate the properties of the half-filling point in lattice QCD (LQCD), in particular the dis¬ 
appearance of the sign problem and the emergence of an apparent particle-hole symmetry, and try to 
understand where these properties come from by studying the heavy-dense fermion determinant and the 
corresponding strong-coupling partition function (which can be integrated analytically). We then add in 
a first step an effective Polyakov loop gauge action (which reproduces the leading terms in the charac¬ 
ter expansion of the Wilson gauge action) to the heavy-dense partition function and try to analyze how 
some of the properties of the half-filling point change when leaving the strong coupling limit. In a sec¬ 
ond step, we take also the leading nearest-neighbor fermion hopping terms into account (including gauge 
interactions in the fundamental representation) and mention how the method could be improved further 
to incorporate the full set of nearest-neighbor fermion hoppings. Using our mean-field method, we also 
obtain an approximate (/i, T) phase diagram for heavy-dense LQCD at finite inverse gauge coupling /3. 
Finally, we propose a simple criterion to identify the chemical potential beyond which lattice artifacts 
become dominant. 


* E-mail: rindlisbacher @ itp.phys.ethz.ch 
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1 Introduction 


Lattice QCD at non-zero baryon density suffers from the so-called "sign problem" [1]. It also displays a 
new type of lattice artifact: because of the Pauli exclusion principle, the lattice cannot support more than 
one quark per site (and per spin, color and flavor). The main goal of this paper is to characterize the effects 
of this new discretization error, so as to better disentangle genuine physics from lattice artifacts. 

We first briefly recall the foundations of lattice QCD in the presence of a chemical potential /i and some 
well-known features of the sign problem which appears when p ^ 0. We then focus on the half-filling 
point (Sec. 2) and the associated fermion saturation. Using the heavy-dense limit (Sec. 3), we elucidate 
the properties of the half-filling regime and analyze, using a mean-field calculation which takes leading 
nearest-neighbor fermion hoppings and gauge interactions into account (appendix B), how far some of 
these findings generalize to full QCD. Finally, using our mean-field method, we obtain an approximate 
T) phase diagram for heavy-dense LQCD, which can be compared with current state of the art complex 
Langevin investigations [3]. 

1.1 QCD at finite Density 

In this paper, we adopt the Wilson discretization of the Dirac operator. For a single flavor at chemical 
potential // it reads: 


DxIa,yJb — ^xy^IJ^ab 
3 

Sx I a ,y Jb 

+S^_2^yil+-f4)abUlljix-Z)e-f^'j , (1.1) 

T^cla.yjb 

where {x, y) are multi-indices labeling different lattice sites, (/, J), (a, h) are color and Dirac indices re¬ 
spectively and U^jj{x) is the (/, J) component of a link variable that points from site a; to (a; -I- F). The 
hopping parameter k is related to the bare lattice quark mass m by k = {2 m + 8)”^, and the boundary 
conditions are assumed to be periodic in spatial directions and anti-periodic in the time direction. 

Using matrix notation, we write (1.1) simply as D{y,), keeping its dependency on the gauge held implicit. 
From the dehnition (1.1), it follows that 


D{fj.) = -f5D\-y*)^5 , 

where y* is the complex conjugate of /i, which shows that 

Det(£i(p)) = J 


( 1 . 2 ) 


(1.3) 


is real if Re(/r) = 0, as 

BetiDiy,)) = Det(75L>l(-/r*)75) = Det(i:)l(-/i*)) 

= Det(D(-p*))* Det(£>(/!))*. (1.4) 

Therefore, if Re(p) = 0, the full fermion determinant for two mass degenerate Havors is manifestly non¬ 
negative, 

D^y)) = Det(i^(p))Det(f2(p)) 

= I)et{D{fj,))I)et{D{ti))* = I)et{D{y) D^p)) , (1.5) 
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( 1 . 6 ) 


and gives rise to a well defined effective fermion action 

S'pff[U] = -Tr{\og{D{^i)D\^i))) , 
which can be used, together with Wilson’s lattice gauge action, 

= U^i^)UAx + n)Ul{x+1})Ul{x) , (1.7) 

^ X 

to sample gauge configurations by Monte Carlo. This generalizes to the case of Nf flavors as long as the 
degeneracy factor for each quark mass is even, i.e. as long as there is always an even number of flavors that 
have the same quark mass. 


1.2 The Sign Problem 


For Re(p) ^ 0, the fermion determinant (1.3) becomes in general complex and loses its probabilistic 
interpretation. In order to continue sampling the gauge field by Monte Carlo, one is then forced to invoke 
reweighting techniques. For example, if the theory contains an even number Nf of mass degenerate flavors, 
one can use |Det(£)(/r))|^''^ = Det{D{fi)) Det{D{—fi)) instead of Det(I?(/r))^^ to sample the 
gauge field and treat the complex phase. 


R[;7] = 


|Det(i7(p))r^ ’ 


(1.8) 


as part of the observables. This is called phase quenching. For Nf = 2, the phase quenched ensemble just 
corresponds to LQCD with an isospin chemical potential. The expectation value of an observable 0\U] in 
the original system can then be obtained as follows: 


^ IV[U]0[U] 

JV[U] 

^ imjOm R[(^] |Det(77(/i))r^ ^ {O R), 

/ V[U] R[17] |Det(i7(p))|^^' (R), ’ 

where {O')^ is the expectation value of 0'[U] with respect to the phase quenched system, i.e. 

JV[U] \Bet{D{p))\^f 

Although the reweighting procedure (1.9) is in principle exact, its applicability is limited by the statis¬ 
tical fluctuations of the phase R[U]. To see this, consider the expectation value of the phase, 

(R),= f = (Lll) 

where A/ = / — /g is the difference in the free energy densities of the unquenched and the quenched 
theory and L and Nt are respectively the spatial and temporal system sizes. For sufficiently large systems, 
A/ becomes approximately independent of the system size and (1.11) shows that (R)^ then decays expo¬ 
nentially with increasing system size. As per definition the modulus of R[(7] remains equal to 1, this shows 
that the fluctuations in the reweighted observable, R[t^]/(R)^, grow exponentially with the system size 
and so does the statistical error. The dependency of the latter on the number of measurements is known to 
be proportional to i/%/# meas.. This means that the statistics required to obtain equally accurate results for 
different system sizes must scale like 

required statistics oc . (1.12) 

The exponential growth of the required statistics as a function of the system size in (1.12) is a manifestation 
of the so-called sign problem which is generic in the numerical treatment of strongly correlated fermion 
systems at finite density. 
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2 Half-filling and Saturation 


As mentioned above, the sign problem is generic in the numerical study of strongly correlated fermion sys¬ 
tems at finite density and appears therefore also in solid state physics when working with e.g. the fermionic 
Hubbard model. There, finite density simulations are often carried out at the so-called half-filling point, 
i.e. at the value of the chemical potential where the fermion density assumes half its maximum value. Due 
to the particle - hole symmetry that exists at this point, the sign problem disappears and it is possible to 
simulate much larger systems than at generic non-zero values of the chemical potential. 

The half-filling point also exists in lattice QCD (LQCD). The reason for this is of course that the quarks, 
as they are fermions, are subject to the Pauli exclusion principle, which implies that on each site, each of 
the 2 X Nf X Nc fermionic states can be occupied only once. For two flavor QCD (Nf = 2, = 3), 

the maximal number of fermions per site is 12 and the half-filling point should therefore correspond to a 
fermion number density of 6. Figure 1 illustrates this for the case of a tho flavor system with an isospin 
chemical potential (i.e. a phase quenched two flavor system): it shows the isospin number density 

niifi) = = J V[U] Det{D{n)) Det(f?(-/i)) (2.1) 

as a function of the chemical potential ^ for different system sizes. The figure also illustrates that the 
location of the half-filling point is independent of the system size as the different curves (corresponding to 
different system sizes) all cross in precisely the same point at n/ = 6, « 1.3. In figure 2 we also show 

the condensates tp') and ('!/’74't3'0) as a function of isospin-/i for the same set of system sizes. Also 
here we see volume independency at half-filling as the different curves cross in a single point at /r « 1.3, 
but in contrast to the situation with the isospin density, the values of the condensates at half-filling are not 
half-way between their corresponding values at ^ = 0 and at saturation /i —oo. 



/J 

Figure 1: Average isospin number density as a function of isospin p for different system sizes, recorded at k = 0.15 
and (5 = 5.0. All curves cross in a single point, the half-filling point at p « 1.3, where the isospin number density 
assumes half its maximum value. 

To see where the volume independency comes from, we define the local isospin density, 


nj{x,p) = 2 Ty.e{iYc 4 ^D{x,x + 4:,p) D~^{x + 

- D{x+\x,p) D~^{x,x-\-2,p)y^ , (2.2) 
where is the trace with respect to color and Dirac indices, such that 


, , 1 dlogZgip.) / 1 ^ ^\ 

= L3]V, 3p = <»'•>')) • 

\ ^ ' q 


(2.3) 
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Figure 2: Quark (left) and "isospin" (right) condensate as a function of isospin fj, for different system sizes, recorded 
at K = 0.15 and j3 = 5.0. All curves cross in a single point, the half-filling point at /r ~ 1.3, but in contrast to the 
situation with the isospin density in Fig. 1, the values of the condensates at half-filling are not in the middle between 
their values at p = 0 and p —^ oo. 


and consider its average space-time variance, given by 

((’*' - ^ 

This quantity is shown on the left-hand side of figure 3 for the same ensembles that were 
As can be seen, (2.4) is essentially zero at the half-filling point, which means that the half-filling is realized 
very homogeneously. As, according to the right-hand side of figure 3, also the ensemble variance 

N, {{n])^ - {m)]) (2.5) 

of the total isospin density (2.3) is approximately zero at the half-filling point, we can conclude that the 
homogeneity of the local isospin density (2.2) at this point is due to the fact that its dependency on the 
gauge field is highly suppressed at half-filling, which explains the volume independency. 


(2.4) 

used in figure 1 . 




Figure 3: Average isospin density space-time variance eq. (2.4) (left) and isospin density ensemble variance (right) as 
functions of p for different system sizes, recorded at k = 0.15 and j3 = 5.0. The figures show that at the half-filling 
point at /i « 1.3, both the space-time variance and also the ensemble variance of the isospin density are essentially 
zero, which suggests that the isospin density becomes independent of the gauge fields for this particular value of p. 

In figure 4 we also show the average sign (or average phase) (1.11) for a two-flavor system as a function 
of the chemical potential (left) and as a function of the isospin number density (right), for two different sets 
of simulation parameters: once for a system with Nt = 4"*, inverse coupling /3 = 5.0 and hopping 
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parameter K = 0.15 and once for one with (5 = 5.3 and/c = 0.175. As can be seen, for 

both systems the average sign has a local maximum at the half-filling point, where it almost reaches a value 
of one again. However, for the larger system a clear deviation from one is visible. The right-hand part of 
Fig. 4 also shows an approximate symmetry about the half-filling point. 


1.0 

0.8 

A 

;§■ 

V 0.4 

0.2 

0.0 


Figure 4: Average sign as a function of p (left) and as a function of the isospin number density (right) for two different 
sets of hopping parameter k, inverse coupling j3 and system size. The figures show that the average sign has a maximum 
at half-filling where its value is close to one. For the larger system there is however a clearly visible deviation from 
unity. 

Although the average sign is not exactly one at half-filling, the fact that it has a local maximum, implies 
that isospin number density and quark number density must be the same at this point, as we have 

K)(/i) - {ni){n) = 

which vanishes at an extremum of (R)^(/r). This is of course also true at the minima of (R)^(Ai), which 
means, that although the sign problem is worst at these points, the quark number density could there be 
obtained with high accuracy. The problem is, however, to determine the location of these points, as they 
will not be volume independent as is the case for the half-filling point. In figure 5, we show (2.6) for the 
same data that was used in figure 4. The density difference was obtained by reweighting the quark number 
density from the quenched to the un-quenched system and then subtracting the isospin number density. 

We mentioned that, in the Hubbard model, the sign problem disappears at half-filling due to a particle- 
hole symmetry. We could now ask if also in LQCD such a particle-hole symmetry is the origin of the 
(almost) disappearing sign problem at half-filling. A first indication that the answer to this question could 
be yes, comes from Fig. 6, which shows, again for a two-flavor system, coupled to an isospin chemical 
potential /r, the quantity nj^sat — {ni){l^) for two different ensembles, where nj^sat = (ni) (n) is 

the saturation value of the average isospin density {nj){iJ.) at large p. We could think of nj^sat — (^/)(/^) 
as being the density of "holes in the saturated state". 

For the isospin density itself we would expect (n/)(/i) oc. ^ for jl, where = — log(2«;) is the 
half-filling value of the chemical potential and Q = 2 is the charge of an elementary excitation (charged 
meson). The fits in Fig. 6 now show that nj^sat — behaves in precisely the same way for ^ fl, but 

with Q = —2, which means that the "hole-excitations" from the saturating state at large chemical potential 
behave just like the anti-particles of the elementary excitations of the theory at small chemical potential. 
The strong interaction treats therefore fermion holes in the saturated state in a similar way as the fermions 
themselves in the vacuum, such that the holes also have to form gauge invariant combinations of the form 
of mesons and baryons. 

Having this symmetry between particles at zero density and holes at saturation, makes it plausible that 
at half-filling, where particles and holes appear with equal densities, they indeed behave symmetrically. 
The same should be true for the non-phase-quenched system where the elementary excitations would be 
baryons with Q = 3 at zero density and holes of the form of anti-baryons with Q — —3 at saturation. 


1 /dlogZ(fj,) dlogZq{fj.)\ 1 d\og{R)^{fi) 

ap dfi ) ~ WWt dfj, 
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Figure 5: Difference between average quark number den¬ 
sity and average isospin number density as a function of /r 
for the same ensembles used in Fig. 4. The quark number 
density was obtained by reweighting. 



Figure 6: The figure shows a log-plot of ni^sat — {ni){fi) 
for the two ensembles used in Fig. 4. The dotted lines cor¬ 
respond to fits of the form n/,sat—(n/)(/r) = 

to the data for > 1.5 p,, where p = — log(2«:) is the 
half-filling value of p in each ensemble. As can be seen, 
both fits yield Q « —2, although the simulation parame¬ 
ters are quite different. 


Again this suggests that at half-hlling, where baryons and baryon-holes have equal density, we should have 
particle-hole symmetry. 

2.1 Dominance of Lattice Artifacts: a Simple Criterion 

In contrast to the systems studied in solid state physics, where the lattice is physical and half-hlling describes 
a physical state as in the fermionic Hubbard model, the lattice in LQCD is jut a regulator which has to be 
removed (e.g. in numerical works, by doing continuum extrapolation) in order to extract real physics. 
The half-hlling state of a LQCD system is therefore just an unphysical lattice artifact and should be of no 
relevance for the continuum physics. 

As long as we are running simulations sufficiently far away from the continuum limit (such that the maximal 
fermion density on the lattice (in physical units) is far below the density where the Pauli exclusion principle 
would become important in the continuum), the hrst appearance of a minimum in the average sign as a 
function of increasing chemical potential (or as a function of increasing isospin number density) as shown 
in hgure 4, indicates the point where lattice artifacts start to become dominant. This becomes clear by 
noting that the average sign is a measure for the overlap between a system with a quark chemical potential 
and the corresponding phase-quenched system (with an isospin instead of a quark chemical potential). As 
the right-hand side of hgure 4 shows, the average sign drops dramatically as soon as the isospin density 
becomes non-zero, rehecting the fact that a pion condensate develops in the isospin system, leading to a 
ground state which is rather different from that of the corresponding system with a quark chemical potential. 
The fact that the overlap between the two systems starts to increase again with increasing p (i.e. that the 
average sign starts to increase again) indicates, that the inhuence of the Pauli exclusion principle starts to 
become dominant on the lattice: the fact that the quark chemical potential favors the u and d quarks, while 
the isospin chemical potential favors u and d, starts to become less important as for example the "d-hole" in 
the system with the quark chemical potential starts to play the role of the d in the isospin system and vice 
versa. This happens at densities signihcantly below half-hlling, and measurements of observables taken at 
larger densities should be considered with caution. 


3 Heavy-Dense Analysis 

In the following section we will analyze the half-hlling point of LQCD in the heavy-dense limit and at strong 
coupling. These simplihcations will allow us to understand the origin of the properties of the half-hlling 
point. 
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3.1 The Heavy-Dense Approximation 

The goal of this section is to obtain an expression for the Partition function in the so-called heavy-dense 
approximation at strong coupling, which can be integrated analytically. The heavy-dense approximation 
usually consists of taking the limit k —> 0, p, —?► oo, keeping n finite while terms proportional to just k 
or Ke~^ vanish. In our case, it would be more appropriate to call it the static quark limit, as we will just 
drop the spatial hopping terms from the Dirac operator (1.1), but not the one corresponding to backward 
hopping in time, which is proportional to Ke~^. The resulting fermion determinant is well-known as the 
leading term in the spatial hopping expansion that was already used in [6] to derive an effective model for 
QCD or in [2] to obtain a holomorphic action for complex Langevin simulations of LQCD. 

Let us start by expanding the determinant of the Dirac operator (1.1) in terms of closed fermion loops 
Cio of length/o [4]: 

Det(D) = exp(Tr(lnD)) = exp(Tr(ln(l — k i/))) 

= expf-Hy Tr(i/'=)J 

^ k^l ^ 

= exp(-^ H 

^ L {c,4"=i 

= expf"^ ^ ti'c,d(ln(l - 
^ {C<„} 

= ll ll det,4l-K‘°McJ , (3.1) 

'0 {C,o} 

with Me, ^ 2,2 ••• ^ 2 ,,^ being the matrix product of hopping terms = (S'-I- along the 

contour Ci^ of length Iq, where S and T are the spatial and temporal hopping matrices defined in (1.1). The 
subscripts "c" and "d" of "det" and "tr" indicate that the operators act on color and Dirac space only. To 
get from the second to the third line in (3.1), we have used that k can be written as k = Iq n, where Iq is 
the length of the contour Cig, and n the number of windings around Cig. The change in the denominator 
(k ^ n instead of fc — >■ Zq n) comes from the fact that Tr(iJ*‘’") contains tTc,d{{Mci^) ) exactly Iq times, 
once for every possible base point (or starting point) of a loop on Cig, which cancels the Iq in the denomi¬ 
nator. 

As already mentioned, the static quark approximation now consists of dropping the spatial hopping 
matrix S from (1.1) and (3.1), which ultimately means that Mc,^ is non-zero only if C/p wraps around the 
temporal direction, i.e. Iq = rit (we use from now on a lower-case rit to represent the temporal system size) 
and has to be one of the ■ riy ■ Uz matrices 

- P, 0 (1 - 74)"‘ = (3-2) 

or their time-reversed versions 

- 0 (1 + 74)”* = (3.3) 

where P^ is the Polyakov loop at spatial position x. The minus signs on the left-hand sides of (3.2) and 
(3.3) come from the anti-periodic boundary conditions in temporal direction. With this approximation, the 
determinant (3.1) reduces to a product of single site determinants: 

Det(P) = ]^(^detc,d(l + «”* Pj 0 (1 - 74 )”*) detc,d(l -I- k”* Pj 0 (1 -t- 74 )”*)) 

X 

= ]j(detc,d(l + (2Ke'")”*P5 (g) ((5 1 ° )) detc,d(l + (2Ke“^)”*P| 0 (^0= °))^ 

X 

= n(det^(l + (2«;e'^)"‘P5) det2(l + (2Ke-'^)”*Pj)) , (3.4) 
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More precisely, the last line of (3.4) shows that for each site, we get four factors, each in the form of a 
determinant in color space; the first two correspond to the two distinct states of a spin-1/2 particle while 
the second two are due to its anti-particle. These color-determinants can be re-written in terms of traces of 
powers of the Polyakov loop; 


det,(l + (2Ke")"‘fl,) = 

■ (» + (1 + + (2Ke'')”*P5)„,, 

3 

= 1 + • • • P-.,a,b, 

A.-1 

= 1 + (2^e'^)-tr,(P,) + (trg(P.) - tr,(P|)) 

+ - 3 tr,(P|) tr,(Ps) + 2 tr,(P|)). (3.5) 

By using that the characteristic equation for a matrix X G SU (3) reads 

X^(A) = detc(Al-X) = A^-trc(X)A^+trc(x1')A-1 = 0 (3.6) 

and, according to the Cayley-Hamilton theorem, X satisfies its own characteristic equation, Xx (^) — ^ 
(in the spectral sense), it follows that 



X^ =1 + %T^(X)X'^ - trc(2f1’)X, 

(3.7) 


= Xl' +trc(X)X - trc(Xl') 1, 

(3.8) 

and therefore 

tre(X^) =3 + trc(X)trc(X^) - tre(X'l') trc(X) , 

(3.9) 


tr,(Xt) =(tr2(X) - tr,(x2))/2, 

(3.10) 

which can be used to simplify (3.5) further; 


detc(l + (2k 

TP.) = 1 + (2Ke'")”‘tr,(P5) + (2Ke'^)""*tr,(P-^) + (2«e'^)^"‘. 

(3.11) 

Finally, expressing the traces in (3.11) in terms of the eigenvalues of the Polyakov loop Pj, 



Ai(s) = A2(s) = A3(3;) = l/(Ai(s)A2(s)) 

(3.12) 


we arrive at 


det^(l + (2«;e^)"*Ps) = -p(2«;e'^)”*)^ 

. +(2Ke^)”‘)"(e-*®==(") +(2Ke'^)"*)" , (3.13) 


and similarly, 

det2(l + (2«;e-'^)"‘P|) = +(2K;e-'^)”‘)^ 

. (p®i(s) +(2Ke-^)”‘)^ (3.14) 

The single site fermion determinants in (3.4) reduce therefore to the simple form; 

Det(D(0i,02;M,K,nt)) = 

(^e‘(^i+ff2) +(2Ke^p)'(e-‘^^ +(2Ke'^)”*)"(e-‘®^ +(2«;e'^)”‘)" 

• +(2Ke-'")”‘)^(e*®i +(2Ke-'")”*)^(e‘®" +(2Ke-^)"‘)^ (3.15) 
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Using the Haar measure for the trace of the Polyakov loop, 


dP = il(6li,6»2)d6iid6»2 , 


(3.16) 


where 


P(01,02) 


24 tt' 




-A, 


i<j 


(sin(0i - 62 ) - sin(20i + 62 ) + sin(di + 202 ))^ 

6 TT^ 


(3.17) 


and neglecting the plaquette action (i.e. considering the strong coupling limit /3 = 0), we can write down a 
heavy-dense partition function for QCD with two flavors, u, d: 


j , TT-i , A) — j dPx Det (Pa; (p-ii ; )) Det (Z^a; (pc/, Tlf )) 

/ f f \ 1^1 

/ / d0id6'2i7(dr,02)Det(P(6ii,6»2;p„,K„,nt))Det(P(6*i,02;Pd, Kd,nt)) j , (3.18) 


.Zsi/iu./Jd.KujKd.flt) 


where K„/d and Hu/d the hopping parameter and the chemical potential corresponding to the u/d flavor, 
A is the spatial lattice with |A| sites, and Zs{idu, Pd, k«, Kd, nt) is the single site partition function. 

The integral (3.18) can be evaluated explicitly but results in a rather lengthy expression. For the purpose of 
illustration its explicit form is shown in appendix A for the case of degenerate quark masses, Ku = Kd = k, 
once with an isospin (p = Pu = — Pd) and once with a quark (n = fiu = Pd) chemical potential. 


3.1.1 Heavy-Dense Staggered Fermions 

The Staggered fermion operator is given by 


ddstx,I;y,J — mj dxydjj -f ‘2jji ^ (P/i P)^a:—p,y) 


y.= l 


(3.19) 


with the Staggered phase Pfiix) = ^_j^^3a:i-i-2x2+a,3 lattice fermion mass m. The full static 

quark fermion determinant for (3.19) is more difficult to obtain than the Wilson fermion analogue (3.4). 
The reason is that Staggered fermions allow for retracting fermion loops, such that even without spatial 
hoppings, several distinct temporal loops are possible. However, if we are just interested in the heavy-dense 
limit, where fi and m are assumed to be large while is of order one, the leading term reads 


Bet{D,t) cx J|dete(l -f P{x)) . (3.20) 

X 


Comparing (3.20) with (3.4) reveals that the Staggered heavy-dense determinant (3.20) can be obtained 
from the heavy-dense limit of (3.4) by replacing 2 k by 1/(2 to) and dropping the square of the color- 
determinant. 


3.2 Analytic Expressions for some Observables 

In this section we derive some observables with respect to the system described by (3.18). 


3.2.1 Average Sign 

In terms of (3.18), with = Kd = k, the average sign is computed as 


{e^'^'*‘){p.,K,nt,nxnynfj = 


Zs{p,,fi,K,nt) 


(3.21) 


^Zsip,,-p,,K,nt) ^ 

where f is the complex phase of the single flavor determinant (3.4), and Ux riy riz is the spatial system size. 
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3.2.2 Isospin and Quark Number Density 

The isospin density is obtained from (3.18) by 




and the quark number density, analogously by 


1 dlogZs{n,-fi,K,nt) 
nt dfi 


j{n,K,nt) = 


1 dlog 


nt 




For the difference rig(/i, K, rii) — n 7 (p,, k, n*), we find 


nq{^,K,nt) - ni{^j,,K,nt) = 


_l _ d\og{e^'^){fi,K,nt,n^ Uy 

Tlx Uy Hz nt dfi 


(3.22) 


(3.23) 


(3.24) 


as in full QCD. 


3.2.3 Heavy-Dense Quark Propagator 

The heavy-dense quark propagator is most easily obtained by using that [ 6 ] 

(1 -KT{fi)y^ = (1 - i^T+{yy^ + (i - KT-{yy^ - i, ( 3 . 25 ) 

where 

Tyy) = U4 {x) 0 (1 - 74 ) e^ (3.26) 

Tyy) = (- 1 )'“'‘’ 4 - 4 ,, Ul (cc - 4) 0 (1 + 74 ) e-'^ (3.27) 

and 

Tx,y = Ty + T^ y . (3.28) 

We then obtain: 


\^._0 ' x,a,I\y,b^J 


(1 - 74)a,b 


00 V 

^5(a^4,y4) Xl(“(2Ke^)”‘P5(2/4))''] + 

fc=0 ' 


( 1 + 74 ). 5 




5-x-y (j + {2 ne^y Py{y4)y^)^ ^ 

, (l+74)a,b^ ^ 

“1“ ^x^y J 5 (3.29) 


where Py(jjA) is the Polyakov loop with base point y and Py{x4, 7/4) is a Polyakov line at spatial position y 
that, for X 4 ^ 2/4, starts and ends at euclidean times 0:4,7/4 respectively (by going in positive time direction), 
while Py{x 4 ,X 4 ) = 1. The inverse of the matrix (1 + (2 k e^)"‘Pg( 7 / 4 )) can easily be carried out (e.g. 
with Cayley-Hamilton and Leverrier-Faddeev), provided the matrix is not singular: 


{l + {2nen"*Py{y4)y^ 


1 

detc(l + {2Ket^yPy{y4)) 
— {2 k e^) 


( 1(1 -f ( 2 k e^y trc(Pg(7/4))) 

"‘P^(7/4) + {2Key^^Plyy4)) , 


(3.30) 
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where the determinant in the denominator is given by (3.11). Similarly we get; 


/ oo X 

\f^—Q / x^a,I',y,b,J 


^(94.-4) f;(-{2Ke-»)"'P|(„))‘) + 

^ ^ k=o ^ 


(1 - 74)a,b 




4.5 ■- * (P|(i/..»4) (1 + (2Ke-»)-pt(i,4))-‘) 

, (l-74)a,6 


I,J 


-S.,ySlJ, (3.31) 


and therefore with (3.25): 

(1 — ^/.y 5 j 

= (1 - 74 ),., {Py{x,,y,) (l + (2«e-)"‘P,(z/4))-\,, 

+ (1+74),., {Pliy 4 ,Xi) (1 + ( 2 K:e-'^)”‘P-^(t/ 4 ))”\,^) • 

(3.32) 


The heavy-dense quark propagator (3.32) could be used to get analytic expressions for average meson and 
baryon propagators from which one could determine the corresponding masses. But due to the projection 
properties of (3.32) in Dirac space, the mass spectrum will be highly degenerate. We will therefore use 
(3.32) in the following just to facilitate the derivation of some other observables. 


3.2.4 Local Isospin Number Density 

The local isospin number density in the heavy-dense case can be obtained, by starting from equation (2.2) 
and replacing the full Dirac operators and propagators by their corresponding heavy-dense counterparts, i.e. 
just dropping the spatial hopping terms from (1.1) and using (3.32) as the corresponding propagator. We 
then obtain: 

ni{x,fi) =2 Re(trc,d|£*(a:, a;-I-4, p) -I-4, X,/r) 

- D{x + 4, X, y.) D~^{x, x -f 4, 

= Re(trc,d|( 2 «:e'")”* {Px{xi) (l -f ( 2 Ke'")”‘P 5 ;(x 4 )) (g) (1 - 74 ) 

- (2Ke“'")"* {pI{xa) (1 + (2/ce“^)"'*P^(x4)) (g (1 -f 74)|) 

= 4 Re(trc|(2Ke'^)”* (P5;(x4) (l -f (2 k e^)"*P5;(x4)) 

- (2^6-'^)”* (Pj(x 4 ) (1 + (2/re-'^)”‘P-^(x4))”^)}) , (3.33) 

where after the second equality sign it was used that (1 ± 74 ) are orthogonal projectors and that for example 

[-ne^ Ua{x) (g (1 - 74 )) • (-(2 Px{xa -f 1 ,X 4 ) g (1 - 74 )) 

= ((2Ke^)”‘Pj(x4)g(l-74)) , (3.34) 

i.e. the link variable Ua{x) coming from the term in the Dirac operator that is not projected out by the 
(1 — 74 ) of the propagator, closes the Polyakov line Px{xa + 1 ,X 4 ) to a Polyakov loop Px{xa), starting 
and ending at x. Analogously: 

{-Ke~^ P|(x-4) g (1 + 74 )) • (-(2Ke“^)"‘ ^pi(x 4 ,X 4 -f 1) g (1 + 74 )) 

= ((2Ate-^)"‘P-^(x4)g(l + 74)) . (3.35) 
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3.2.5 Condensates 

The quark condensate for a two-flavor system, 

= (uu) + (dd) , (3.36) 

is in the static-quark limit obtained by taking the expectation value of 

— (trc.d.t((l - KuT(/i„))_^) -|-trc,d,t((l - KdT{fid))^^)) 

Tit ^ ' 

= ^(trc,d((l + (2K„e'"")"*Ps) ^0(1-74)) -f trc,d((l + ^0(1-1-74))) 

+ + ( 2 Kde''‘*)"‘Pa;) ^ 0 ( 1 - 74 )) + trc,d((l + (2Kde“'"‘^)"*P^) ^0(1+74))) 

= 2 (tr,((l + ( 2 /c,e^“)”*P 5 )”^) + tr,((l + (2 

+ 2 (tr,((l + (2Kde'^‘*)"‘Ps)”^) + tr,((l + (2, (3.37) 

where (l — kT(p))_ ^ is the static-quark propagator (3.32) for y = a; and Kd and p„, Hd are the hop¬ 
ping parameters and chemical potentials for the two flavors. 

Analogously the isospin condensate, 

(^74T3'(/)) = (u74u) - (d74d) , (3.38) 

is obtained as the expectation value of 

— (trc,d,t((lc.i ® 74)(l - KuT{fiu))^^) +trc,d.t((lc,i ®74)(l - Kd P(/4d)) ^ 

= 2 (^'^c,d((l + (2 Ku e^“) P 5 ) 0 (74 ~ 1)) + trc,d((l + (2 e P^) 0(74 + !))^ 

— ^(trc,d((l + (2Kde'''*)"‘P5;)”^ ® (74 - 1)) + trc,d((l -f (2 ® (74 + 1))) 

= 2 (-tr,((l +( 2 K„e'^")”‘P 5 )-^) + tre((l + (2 e-'^“)”‘Pj)-')) 

- 2 (-tr,((l + (2Kde'^'‘)”‘Ps)”') + tr,((l + (2e-^'‘)”‘pi)”^)) , (3.39) 
where in the first line is the identity with respect to color and time indices. 

3.3 Half-filling Point 

Having obtained the heavy-dense partition function for two-flavor LQCD in the strong coupling limit, and 
defined some basic observables, we now use these results to study the origin of the properties of the half- 
filling point in the heavy-dense/strong coupling limit, and check, how far our findings could also apply to 
full LQCD. 
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3.3.1 Location of the Half-filling Point 

To find the location of the half-filling point, we consider the isospin density (3.22) in its integral form: 


= — ^log( [ dP\I)et{D{fj.,K,nt))\' 
n-t \J 


1 




1 


4 / 


dP 


dP 


d\Bet{D{n,K,nt))\ 

dfi 

a|det^(l -h (2«;e'^)"*P) det^(l -h (2Ke-'")"‘P^) | 


dfi 


\ -i\ 


, — / dP Re((2Ke'")”*tre(P(l + (2Ke^)"‘P)’ , 

Zs[fi,-fi,K,nt) J 1 

-(2/te-'^)"‘ tr,(P^ (1 + (2Ke-^)"‘P^)“^)} 

X |det^(l + ( 2 Ke'")”‘P)det^(l + ( 2 «:e-'^)”*P^)|^ ( 3 . 40 ) 


For sufficiently large values of fi, the second term within the curly brackets in (3.40) can be neglected; 
setting fi = — log(2K), this term is of order C2((2 k)^"*) while the first term is of order 1 and becomes 
completely independent of P, as 

Re{tr,(P(l+P)“')} = l{tr,(P(l+P)“') + tr,(P^ (l + P^)”^)} 

= ^{tr,(P(l+P)-') +tr,((P+ 1)-')} 

= ^tr,((l+P)(l-hP)-') = ^, (3.41) 

and we obtain for the isospin density (3.40): 

n/(/r =-log(2K),K,nt) = 6 + C>((2 k)^”‘) , (3.42) 

which means that fi = p. = — log(2K) corresponds to the half-filling point, independently of rit (up to 
Cl((2/c)"‘) corrections), as is visualized in figure 7, where the isospin number density (3.22) is shown as a 
function of fi for different system sizes and it can be seen that the curves, corresponding to different system 
sizes, all cross at the half-filling point, just as in the case of full LQCD shown above in figure 1. From 
(3.40) and (3.41) it also follows that the dependency of the local isospin density (3.33) on the gauge field 
is highly suppressed at half-filling: gauge field dependent terms are suppressed by a factor of C)((2«;)^"‘) 
such that the space time variance of (3.33) nearly vanishes at half-filling, again just as in the case of full 
LQCD shown above in figure 3. 


Nevertheless, in full LQCD, the half-filling point is shifted towards a slightly larger value of p,: while 
in the heavy-dense case, we have p = — log(2 k) « 1.2 for k = 0.15, we find in full LQCD, for the same 
value of K and with /3 = 5.0, that the half-filling point is located at /i « 1.3. This shift is due to spatial 
fermion hopping which causes corrections to the quark mass. Corrections coming from spatial hoppings and 
a finite value of /3 will of course also appear in the expressions for observables and give rise to deviations 
that could be much larger than C2((2 These corrections are considered in Secs. 3.5 and 3.6. 
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Figure 7: Isospin density in the heavy-dense approxi¬ 
mation as a function of ^ for nt € {4, 8,16,32} and 
K = 0.15. Compare with Fig. 1. 



Figure 8: Isospin density space-time variance in the 
heavy-dense approximation as a function of ^ for rit € 
{4, 8,16, 32} and k = 0.15. Compare with Fig. 3. 


3.3.2 Particle-Hole Symmetry and Average Sign 

Consider the single flavor fermion determinant (3.15) under the transformation /r —> 2 jl — ^ that corre¬ 
sponds to reflection about the half-filling point: 

Det(Z7(2/x —/r, K, nt)) = detK^l + {2Ke^)~^*P) det^(l + ((2 k)^ e'^)"'*P^) 

= (2Ke^)”®”‘det^(P)det^(l + (2ne^)”‘P^) det^(l-F ((2 k)^ e'^)"*p1') 

= (2«;e'')”®"‘ det^(l + (2«:e'")”*p1') (l -F C>((2 k)”*)) , (3.43) 
where the last line of (3.43) is true for fi G [0, 2 fi]. Comparing (3.43) with 

Det(Z7(^, tv, nt))* = Det(Z7(—/i, k, nt)) = det^ (l + (2 k e“^)"‘P) detc(l + (2Ke^)"*P^) 

= det^(l + (2Ke'")”*p1') (l + C>((2 k)"*)) , (3.44) 

we see that they agree up to an irrelevant gauge field independent pre-factor (which is unity at the half¬ 
filling point) and terms which are suppressed by at least a factor of (2k)"*. As the complex conjugate of 
the fermion determinant can be interpreted as resulting from a charge conjugation of the original fermion 
fields, this shows that the half-filling state possesses an exact T = 0 particle-hole symmetry also in the 
heavy-dense/strong coupling limit of LQCD. 


We can also check the behavior of the heavy-dense quark propagator (3.32), 

= (1 - 74),,, {Py{x,,y,) (l + {2 Py{y,))-^) 

+ (1+74 ),., (pJ( 2/4,X4) (1 + (2Ke-'*)”‘Pj(2/4))”'),,^) , 

(3.45) 


under the transformation /i —)■ 2 ji — p,. For 7/4 yF X 4 and p S [0, 2 p] we find: 

(1 - «T(p));l^,^^_, ^ ^ _(1 + o((2k)^)) . (2Ke^)“°'^(--^-"*)(-l)^(«^>-) 

(1 - 74),., {P^{y4,X4) (1 + (2Ke'*)"*Pj(y4))”')4„4 (3.46) 

and 


(1 


— K 


n-y)) 


-1 

x,a,I;y,b,J 


-(1 + 0((2k)")) • (2Ke'*)”°'^(^^-""’"‘^(-l)"("">""^ 

(1 + 74),., { Py { x 4 , y 4 ) (1 + (2Ke^)"*P^(y4))”')4,j ’ (3-47) 
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which by noting that 


(Py{x4,y4) {1 + i2Ke^)'^*Py{y4)) ^) = ((l + (2 Ke'")”*Py(a:4)) ^ Py{x4,y4)) , (3.48) 

and remembering that 

75 (1 ± 74 ) 75 = (1 =F 74 ) , (3.49) 

shows that: 

(1-«:T(2/i-/i)) = -75 (1 - «:T(-^)) 75 (1+ C>((2«:)^)) 

= -(1-kT(^))^(1 + 0((2k)')) , (3.50) 

and 

{1 - kT{-{ 2^-y))) = -75 (1 - Kr(/r)) 75 (1+ C>((2k)^)) 

= -(1 - nT{-y)y (1 + 0{{2Kf)) , (3.51) 

i.e. up to a minus sign and terms of higher order in {2k), we find that for j/4 ^ X 4 , reflecting y, about y turns 
the quark propagator into its Hermitian conjugate. For the equal time part, 2/4 = X 4 we find in a similar 
way: 

(1 - ^ 5 ((J - 74,,. («w - {I + (iKe-rptM)-]) 

+ (l+74)a^U^i,J + 0((2«)'"‘(2K:e'^)"‘))) , (3-52) 


and 

(1 - j -)■ -(^(1 + 74 ) 0,5 (<^/+ - (1 + (2Ke'")”*Pjj(y4))j j) 

+ (1- 74 ) 0,5 (<5/+ + 0((2K)"”‘(2K:e'^)”‘))) . (3.53) 
For the two condensates introduced in Sec. 3.2.5, we then find for the phase quenched case: 

(^’'0),(2A-/4) = 12 +(^74T3^/i)^(/i) + 0((2 k)"*) , (3.54) 

and 

{^j 4 T 3 'ijj) ^{2 y - y) = {^tlj)^{y) - 12 + C>((2k)”*) , (3.55) 

which is illustrated in Fig. 9. In contrast to the situation in full LQCD shown above in Fig. 2, in the static 
quark/strong coupling limit, the values of the condensates are at half-hlling exactly in the middle between 
their corresponding values at y = 0 and for y ^ 00 . 

In the un-quenched case, the relations would be: 

{^ip){2y - y) = 12 + (^74V')(/r) + C>((2tt;)”*) , (3.56) 

and 

(^74V’)(2/i - ri) = (^’V’)(f) - 12 + C>((2k)”*) . (3.57) 


Returning to the single flavor fermion determinant (3.15), we can see from (3.43) and (3.44) that it is 
essentially real at y = y and the average sign (3.21) should therefore be 1 + C9((2 k)"*). This can also be 
verified more directly by using (3.15) to show that for y = y, we have: 

dete(l + (2Ke'^)"‘P) = 1 + (2«:e^)"‘tre(P) + (2«;e'^)2"‘tre(P^) + (2«:e^)3"‘ 

= 2 + tr,(P) + tr,(pl) = 2(1 + Re(tre(P))) , (3.58) 
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Figure 9: Quark (left) and isospin condensate (right) in the static-quark limit, defined as the expectation values of (3.37) 
and (3.39) respectively, shown for Ku = Hd = 0.15 as functions of an isospin chemical potential jj. — fiu = for 
rit £ {4,8,12,16}. Compare with Fig. 2. 


i.e. we find again that the dominant part of the single flavor fermion determinant (3.5) at half-filling, 

det^(l+P) (1 + C>((2 k)2"*)) (3.59) 

is real and positive up to terms of order 


The average sign (3.21) in the heavy-dense limit is shown in figure 10 as a function of the chemical 
potential (left) and as a function of the isospin number density (right) for the same system sizes and values 
of the hopping parameter k that were used to generate the corresponding plots in full LQCD shown above 
in figure 4. By comparing the two figures, it can be seen, that in the heavy-dense case, the average sign is 
much more symmetric about the half-filling point and at the half-filling point itself, it deviates much less 
from being unity than in full LQCD, which is again due to the absence of spatial fermion hoppings, such 
that the deviation is really just of order C9((2/c)^"‘). 




Figure 10: Average sign (3.21) in the heavy-dense/strong coupling approximation as a function of the chemical potential 
p (left) and as a function of the isospin number density (right) for the same sets of hopping parameter and system size 
that were used in figure 4 above for full LQCD. 


3.4 Mean-Field Method for Heavy-Dense QCD at Finite Gauge Coupling /3 

So far we have only considered heavy dense LQCD in the strong coupling limit, /3 = 0. In this section, we 
would now like to get a glimpse of finite /3 effects by introducing an effective nearest-neighbor Polyakov 
loop action, which is motivated by the leading terms of the character expanded gauge field Boltzmann fac¬ 
tor, as derived for example in [6], and then using a mean-field approximation. Our mean-field treatment is 
similar to one of the approaches considered in [8] but differs from the mean-field calculation in [7] as we 
are considering the Polyakov loops as our effective degrees of freedom while in [”], this role was taken by 
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spatial links. 


The nearest-neighbor Polyakov loop action is obtained as follows: we start with a character expansion 
of the SU(3) Yang-Mills Boltzmann factor (c.f. [5], chapter 3.4), 

- HIE drCr{l3)xr{Up)\ = n{l + E drar{P)Xr{Up)\ , (3.60) 

where the product is over all plaquettes p, Xr{Up) is the character of the SU(3) group element Up (cor¬ 
responding to the plaquette p) in the representation r, and aj.(/3) = Cr{P)/Co{f3) with Co(/3) being the 
expansion coefficient of the trivial character Xo{Up) = 1- We now drop the summation over r and keep just 
the terms corresponding to the fundamental representation. After integrating out the spatial links, we are 
then left with [6]: 

n{l + ^?(/d){xf{P^)XfiPj) + Xf{Phxf(.PJ))} 

{i,3) 

= n{l + a;;‘(/3)(tre(P.) ir,{P]) + ir,{Pl) tr,(P,))} , (3.61) 

{i,j) 


where the products run over all pairs of nearest-neighboring sites, (i, j). If we now require that the Boltz¬ 
mann factor corresponding to the desired effective nearest-neighbor Polyakov loop action should, to lowest 
order, be proportional to (3.61), we find that the effective action should be given by: 

- 5e// = af (13) Y, L*j + Lt Lj) . (3-62) 

(i.i) 


where Li = tidPi)- At this point it should be mentioned that in contrast to the situation with the U(l) 
or SU(2) Yang-Mills action, where the coefficients of the character expansion can be written in closed 
form in terms of Bessel functions, no closed form is known for these coefficients in the case of SU(3). To 
compute af{/3), we made use of the power series representations for Co(/3) and Cf{/3), given in [5] up to 
order i.e.: 


co(/?) 


l + (/3/6)' + 

23(/?/6)^° 

2400 


{/3/6f , (/3/6)^ 


29(/3/6) 

7560 


2 

11 


-f 



2629(/3/6)^^ 

1814400 


n{(3/6f ll(/3/6)" 139(/3/6)" 19(/3/6)^ 

72 120 2880 864 

1241(^/6)^^ 17449(/3/6)^" 

2419200 101606400 ■" ^ ' 


and 


c/(/3) = (/3/6)+(^+(/3/6)^+ 

319(/3/6)^° 2629(/3/6)^^ 

15120 302400 


5(/?/6)^ l3{(3/6f 77{(3/6f 

8 ^ 24 ^ 240 

16133(/3/6)^^ , 17449(^/6)^^ 
4838400 14515200 


-f- 

-f 


139(/3/6)" 

720 

35531(/3/6) 

87091200 


19(/?/6)^ 23(/?/6) 
192 480 

14 

- + .... (3.64) 
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We now proceed by applying the mean-field approximation to (3.62), i.e. we write Li = L + 5Li, keep 
only terms of order 0{5Ld and then write again 6Li = Li — L (and proceed analogously for L*), which, 
after dropping constant terms, leads to the following single site mean-field action: 

- Smf,^ = 6 af (/3) (L* L + L* L^ , (3.65) 

which we use to define an additional probability weight, 

w{L,L*,L,l3,nt) = <^7 ^ 3 ) d " ^ ^ (3.66) 
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that has to be included in the single site partition function defined in (3.18), i.e.; 


2 tz 27r 


Zs{^J,u,^J■d,^^u,Kd,nt,L,L*) = J J d6id92^H{9i,e 2 )Bet{D{ 0 i, 62 -, Hu, Ku,nt)) 

0 0 


j:)et{D{9i,92;Hd,Kd,nt))w{L,L*,L{9i,92),l3,nt) > ■ (3.67) 


For general values of Hu, Hd, Ku and Kd, the product of the two fermion determinants in (3.67) is usually 
complex. As a consequence (L) and {L*) differ and it is a subtle issue how to proceed in a mean-field 
treatment [8]. To bypass the subtleties associated with L* ^ L when setting L = (L) and L* = (L*), we 
define L to be the mean-field value of L with respect to the phase quenched system, where (L*)^ = (L)^ 
holds and the identification L = {L)^, L* = {L*)^ therefore leads to L* = L. The mean-field L is therefore 
determined by solving the following self-consistency equation: 


L = {L)^(hu, Hd, Ku, Kd,nt, L) 

27r 27r 

= ■—- - -— [ [ d9^d92\m,92)H{9u92) 

^00 

\Det{D{9i,92-,Hu,Ku,nt)) T>et{D{9i,92-,Hd,Kd,nt))\w{L,L(9i,92),l3,nt) \ , (3.68) 
where z(;(L, L(0i, @ 2 ),/3, nt) = w[L,L,L{9i, 92), P,nt) and where 

27r 27r 

Zs,q{Hu,Hd,Ku,Kd,nt,L) = J J d9id92l^H{9i,92)\Bet{D{9i, 02; Hu, Ku,nt)) 

Bet{D{9i,92; Hd, Kd,nt))\w{L,L{9i,92), (3,nt) \ (3.69) 


0 0 


is the phase-quenched partition function. After one has determined the stationary L with respect to (3.68) 
for a particular set {hu, Hd, Ku, Kd,nt) of parameters, one can reweight to the non-phase-quenched system 
by using L(^Hu, Hd, Kd, nt) in (3.67) when computing expectation values of observables. 


Of particular interest among the observables that can be measured for the phase-unquenched system, 
are the average Polyakov loop, 

{L)[hu, Hd, Ku, Kd,nt) 


2tt 277 


yy J J d9,d92l^H9„92)H{9,, 02 ) 
0 0 


Zs i^Hu, Hd, Ku, Kd, rti , T (^Hu, Hd, Ku, Kd, Tit 
T)ei{D{0i,02-,Hu,Ku,nt)) Det{D{9i,92-, Hd, Kd,nt)) w{L{hu, Hd, Ku, Kd,nt), L{9i,02), l3,Tit) | , 

(3.70) 


and the average complex conjugate Polyakov loop, 

(T } (^Hu, Hd, Ku, Kd, Tlt^ 

277 277 

1 

(Mu ; l^di 1 ^di 1 ^ (^f^u 1 f^d t 7 ^d7 ) 

Det{D{9i,62] fiu7 f^u7nt)) Bet{D{6i,92] fid7 i^d7nt)) w{L{fj,u7 fJ'd7 ^U7 



0 0 


L%9i,92)H{9i,92) 
Kd,nt),L{0i,92),(3,nt)^ , 

(3.71) 
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for which we now have (Z/ ) ^ (-^) (^f^m ('^) (Mu? ^dj ^t) 

due to the in general complex product of fermion determinants in (3.70) and (3.71). This difference between 
(L*) (p„, p,d, Ku, Kd, nt) and (L) (p„, k„, k^, rij) can clearly be seen in Fig. 11 where we show the 

two quantities for the mass degenerate case (i.e. = Kd) as a function of a quark chemical potential 

/r = = fid- 


average Polyakov loops <L>, <L*> 



Figure 11: Comparison of the average Polyakov loop and average complex conjugate Polyakov loop for a system with 
Tit = 4, Ku = Kd = 0.15 and /I = 5.0 as a function of a quark chemical potential p = = fid, obtained according 

to (3.70) and (3.71). The average Polyakov loop (black curve) initially decreases with increasing chemical potential. 
Such a behavior has also been found in [10] where in their figure 2, the expectation values for Polyakov loop and 
anti-Polyakov loop are shown for a SU(3) matrix model with a (quark) chemical potential. 


It has been mentioned in [8] that in a mean-held treatment, the reweighting method can be no more 
than an approximation scheme: while in a full Monte Carlo simulation, reweighting is exact in the limit of 
inhnitely many sampled conhgurations, a mean-held treatment relies only on the most dominant conhgu- 
ration. In our opinion this argument does not hold as in a mean-held treatment, the expectation value of 
an observable is determined on the active site, usually by integrating exactly over all possible values of its 
conhguration variables, and reweighting should therefore work perfectly well! 

Setting L* = L, the integrals in (3.67) and (3.69) can be carried out exactly, leading to solutions in 
terms of inhnite sums of modihed Bessel functions of the hrst kind; 


fd'd, K-ii, Kd^O/t, L'j — 


/ > fdd, Ku, Kd,nt) 2_^ k\ik-l)\ ^ 

h,l2 k=0 m=0 


2k-h 
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(3.72) 


and 
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k\(k-h)\ ^ 
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I,^k-i2-m){l2afil3) L) 

(3.73) 


respectively, where the coefficients (p„, fid, Ku, Kd, Ut) are given by 


H[9i,e2)T)et{D[9x,e2-, k-u, Ku,nt)) I4ei{D{9i,92-, kd, Kd,nt)) = 

Y, Ci,,i2{ku,kd,Ku,Kd,nt) . (3.74) 

h ,h 
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Unfortunately, the evaluation of (3.72) and (3.73) to the required accuracy is numerically rather expen¬ 
sive and it turned out to be more efficient to use direct numerical integration to solve for L in (3.68) and for 
computing observables. Expressions (3.72) and (3.73) were therefore merely used for cross-check purposes. 

An expression for (3.72) in the more general case of L* ^ L could be found along the lines of [9] where 
the mean-field action (including correction terms) is derived for the case of a pure SU(A^) Polyakov line 
model with a chemical potential. 

3.5 Mean-Field Heavy-Dense Phase Diagram and Finite fi Effects 

In figure 12 we show the average Polyakov loop for the mass degenerate two-flavor case as a function of 
isospin chemical potential p, = /j,„ = —and inverse coupling (3 for five different temperatures, as 
obtained with our mean-field method described in the previous section. As can be seen, for sufficiently large 
f3 a deconflnement transition occurs and the value of f3 at which this transition happens becomes smaller 
with increasing temperature (decreasing rit). 



Figure 12: Average Polyakov loop as a function of isospin chemical potential p = p„ = — p^j and inverse coupling jS 
for five different temporal system sizes: n* = 3 (orange), nt = 4 (blue), nt = 5 (green), nt = 6 (red) and nt = 8 
(purple), computed according to (3.68) with Ku = Kd = 0.15. The two figures show different viewpoints on the same 
data. 

A finite value of the inverse gauge coupling /3 has also an effect on the average sign, which is illustrated 
in figure 13. There we show mean-field results for the average sign as a function of p = p„ = p^ for a 
mass degenerate (k — 0.15) heavy-dense two-flavor system of spatial size U = = 4^ and temporal 

extent Nt = 2, 3, 4, 8. The left-hand part shows the strong-coupling limit, /3 = 0, while the right-hand 
part shows the /3 = 5.0 case. For static quarks, the location of the half-filling point is not affected by the 
finite /3 value, but the overall temperature dependency of the average sign clearly changes: at finite /3 the 
sign-problem becomes weaker with increasing temperature and almost disappears after the deconfinement 
transition. 

Note that the computation of the average sign with our mean-field method comes with a slight complication: 
the active site couples to its six nearest neighbors and there are therefore 6 interaction terms entering the 
mean-field action. But the ratio of links to sites should be 3 : 1 on a periodic 3-dimensional lattice, which 
means that by using simply a formula of the form of (3.21) to compute the average sign with our mean-field 
partition functions (3.67) and (3.69), i.e. 



(3.75) 


we would get a wrong result. To get the correct answer, we have to include another reweighting factor in 
(3.75), coming from one of the passive sites (a pair of one active and one passive site, together with 6 links 
connecting an active and a passive site, describes the smallest representative subset of a larger system, in 
which for example even sites are active and odd sites are passive or vice versa, and the computed sign will 


correspond to such a system). This factor is obtained by taking the ratio of the average fermion determinants 
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in the non-phase-quenched and phase-quenched case. The average sign then becomes: 

Zs{^i, K, nt,L{p, K, nt)) (Det{D) Det(i:>)) (p, k, rit, L{fi, k, rit)) 

Zs,q{iJL, K, nt,L{n, K, nt)) (|Det(i:)) Bet{D)\) k, nt, k, nt)) 

where Det(Z3) is the static quark determinant (3.4) and the factor of 1/2 in the exponent in (3.76) is there 
because the term inside the bracket is now the reweighting factor for two sites. 





Figure 13: Mean-field results for the average sign as a function of p = for a mass degenerate (k = 0.15) 

heavy-dense two-flavor system of spatial volume V = and temporal extents = 2, 3, 4, 8. The left-hand 

part shows the strong-coupling limit, 13 = 0, while the right-hand part shows the (3 = 5.0 case. 

Figure 14 finally shows the average Polyakov loop as a function of temperature and isospin (left) or 
quark (right) chemical potential i.e. it essentially shows the (T, p) phase diagram of our simplified model 
in the phase-quenched and un-quenched case: one can read off how the pseudo-critical temperature for the 
deconfinement transition changes as a function of the chemical potential. The right-hand part of Fig. 14, 
corresponding to the un-quenched case, can be compared, e.g. to the phase diagram in [3, Fig. 2], obtained 
by complex Langevin. In order to simplify the comparison. Fig. 15 shows the data from Fig. 14 using the 
same scales used in [3, Fig. 2], assuming that the lattice spacing of a = 0.15 fm, determined in [3], should 
apply equally well to our system, as we used the same simulation parameters k = 0.04, (3 = 5.8. In our 
effective model, we take the gauge field only in the fundamental representation into account. But at high 
temperature, effects coming from higher representations are much less suppressed than at low temperature 
(the coefficients ar(/3), if included in the effective Polyakov loop action, would appear only with small 
exponents, leading to a weaker suppression of the higher order terms), which explains why in our figure the 
deconfinement transition is generally shifted towards larger temperatures compared to [3, Fig. 2]. 

/S=5.8 /3=5.8 



Figure 14: Average Polyakov loop as a function of temperature 1/Nt and isospin (left) and quark (right) chemical 
potential respectively at /) = 5.8, computed according to (3.68) and (3.70) respectively with Ku = Kd = 0.04. 
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0 . 


Figure 15: Average Polyakov loop as a function of temperature T and normalized quark chemical potential where 
/t = — log 2 k is the half-filling value of in the heavy-dense limit, with k = 0 . 04 . It is the same data as in the 
right-hand part of Fig. 14, expressed in the units used in [3, Fig. 2], in order to allow for a better comparison. 


3.6 From the Heavy-Dense to the Full Fermion Determinant 

The heavy-dense fermion determinant is the leading term in the so-called spatial hopping expansion of the 
full fermion determinant which, for one flavor, can be obtained as follows: 

Det(i:)(/r)) = Det(l - kT - Ks S) = Det(l - k,S{1 - kT)~^^ Det(l - kT) 

-E yTr((5(l-KT)-^)") l.Det(l-Kr), (3.77) 

where S and T are the spatial and temporal hopping terms as defined in (1.1) and we have added a subscript 
s to the spatial hopping parameter, just for book-keeping purposes. The matrix (1 — kT)~^ is the heavy 
dense quark propagator given by (3.32). We can now proceed in a similar way as in (3.1) to write (3.77) in 
terms of a product of smaller determinants of "loop matrices", but this time, the loops are purely spatial and 
the matrices carry also a time-index, i.e. 



Det(i9(p)) = exp< 

^ oo 

-E E E 

tTc,d,t{{McS) 1 • Det(l -kT) 


{Cso}-=l 

= exp< 

[E E 

(log(l - k^oMcso))! • Det(l - kT) 


^ {Cso} 

II 

detc,d.t(l 

•Det(l-KT), 

^ So 


/ 


where the subscripts in trc,d,i, detc,d,t indicate that these operators act on color, Dirac and time indices. 
CgQ describes a closed spatial path of length sq (where, in contrast to (3.1), backtracking is now allowed) 
and the matrix Mcs^ given by the (ordered) product of the matrices 5'(1 — kT)~^ along that spatial path 
(see Fig. 16), 


Me. 


0 


n mod so) (1 '*^)s({i + l)modso) ’ 


(3.79) 


i=0 


where each of the Sxi,xi+i (1 — is understood to be a matrix with color, Dirac and time indices. 
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Figure 16: The figure illustrates the meaning of the spatial loop matrix, introduced in eq. (3.79), for the case of a 
contour Csq with so = 4, which is spanned by the four spatial sites xo, Xi, X 2 and X 3 . In the example, the loop starts 
and ends at the spatial site xg, but as the matrices always appears as detc,d,t(l — ) (or inside a trace 

tfc.d.t {Mc„g )), it is clear that the base point of the spatial loop is arbitrary. 

If we look at a particular matrix element of Me, with time indices X 4 ,,s and x,^^, it corresponds to the sum over all 
possible ways to draw a quark world-line (example depicted in red), connecting (iro, 0:4,5) and (So, XA,e) by winding 
exactly once around the spatial (oriented) loop C 4 (depicted in blue): it starts with a spatial hop from (xo, 0 : 4 , 5 ) to 
(xi, X4,s). The next step involves the static quark time-propagator, which is the sum over all possible ways to go from 
(xi,X4,s) to let us say (xi,X4,ij) (including all possible positive and negative windings around the time direction). 
It continues with another spatial hop from (xi, X4,i^ ) to (x2, X4,i^), followed by another static time-propagation from 
(x 2 , X4,ii ) to (x 2 , X4,i2), and so on. The final step consists of a spatial hop from (xa, X4,i3) to (xo, 2:4,13), followed by 
a final static time-propagation from (xo, 2 : 4 , 13 ) to (xq, X 4 ,e). The full matrix element {Mc,)x, ^ (color and Dirac 

indices suppressed) is obtained by summing over all intermediate time indices X 4 ,ij, X 4 ,i 2 and X 4 ,i 3 . 

Similarly, a matrix element of the n-th power of Me,, e.g. {Me^)x,^s,x,,c^ corresponds to the sum over all possible 
ways to draw a quark world-line, connecting (xo, X 4 , 5 ) and (xg, X 4 ,e) by winding exactly n times around the spatial 
loop (74. 

The final form of (3.78), which could be called a spatial loop expansion, is particularly useful if one is 
aiming towards including spatial fermion hoppings into the mean-field treatment. In this case one can just 
restrict the product over sq in (3.78) to the factors with sq = 0, 2 but still get the full interaction between 
nearest-neighboring sites. In contrast, in the original spatial hopping expansion (3.77), a nearest-neighbor 
truncation would limit the accuracy to 0{k^). 

In order to see how the presence of spatial hoppings changes the properties of the static-quark system, 
it is sufficient to take just the lowest order terms in Ks of (3.78) into account. These terms lead, after 
integrating out the spatial links, to corrections up to order O [kI {(3)) (see appendix B). Using this 
improved effective nearest-neighbor Polyakov loop action within the mean-field framework described above 
in Sec. 3.4, we can for example compute the average Polyakov loop for a mass degenerate two-flavor 
system. The result is shown in the left-hand parts of the figures 17 (phase-quenched case) and 18 (non- 
phase-quenched case), together with the corresponding data obtained from Monte Carlo simulations of full 
LQCD. More precisely, the left-hand part of Fig. 17 shows the full LQCD result together with mean-field 
results obtained with three different effective nearest-neighbor Polyakov loop actions: 

1. the one just mentioned, which takes fermion hopping into account up to order Oi^til Oj* (/3)) (dashed 
black curve), 

2. one that takes fermion hopping into account only at order 0[k^) (dotted black curve), 

3. and the one from Sec. 3.4, which neglects fermion hopping completely (dotted red curve). 

As can be seen, allowing for spatial fermion hopping in the effective theory shifts the value of p, where the 
average Polyakov loop is maximal closer to the corresponding value found in full LQCD. With increasing 
approximation order also the magnitude of the average Polyakov loop gets closer to the LQCD result. In 
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the right-hand part of Fig. 17 we tried to compensate for the truncated spatial nearest-neighbor interaction 
in the effective theory by increasing the value of the inverse gauge coupling /3. Setting /3 = 6.035 in the 
effective theory, the average Polyakov loop obtained by mean-field with the a^*(/3)) effective action 
from appendix B then almost coincides with the corresponding LQCD result for /3 = 5.0. 




Figure 17: The left-hand figure shows the average Polyakov loop for a phase quenched, mass-degenerate two-flavor 
system with UxTiyUzUt = 4^, k = 0.15 and P = 5.0. The solid curve guides the eye through the data points 
coming from a full Monte Carlo simulation, while the dashed line corresponds to the result of a corresponding mean- 
field calculation with our effective nearest-neighbor Polyakov loop action (see appendix B), including terms up to order 
{P)). For comparison, the dotted curves show the corresponding results when using only (7 (k^) (black dotted 
curve) or O (k2 ) (red dotted curve) terms in the fermionic part of the effective action, which results in larger deviations 
from the Monte Carlo data. The right-hand figure shows the same Monte Carlo data (solid curve), but this time it is 
compared to mean-field results obtained at a slightly larger value of the inverse gauge coupling P = 6.035 (trying to 
compensate for the neglected piece of the full gauge interaction). The dashed line, corresponding to the O (P)) 

improved effective nearest-neighbor Polyakov loop action, then almost coincides with the Monte Carlo data. 


1.5 


1.0 


0.5 


0.0 


Figure 18: Average Polyakov loop (black) and complex conjugate Polyakov loop (red) in a non-phase quenched two- 
flavor system with nxnyTiznt = k — 0.15 and P — 5 . 0 . The solid curves are to guide the eye through the data 
points coming from a Monte Carlo simulation of full LQCD, while the dashed lines in the left-hand figure show the 
result of a corresponding mean-field calculation analogous to the one that lead to Fig. 11, but using an effective action 
that takes terms up to order {P)'j into account (see appendix B). For comparison, the right-hand figure shows 

again the data from Fig. 11 (dotted lines) together with the full LQCD Monte Carlo data. 

Finally, figure 19 shows the average sign (3.21) as a function of /i, again for a mass-degenerate (k = 5.0) 
two-flavor system of spatial size V — = 4^, computed by mean-field, using our improved effective ac¬ 

tion and (3.76). The left-hand part shows results for three different values of the inverse gauge coupling 
/3 = 0, 0.5, 0.7 for fixed Nt = 4, while in the right-hand figure, results for different values of Nt = 2, 3,4 
are shown for /3 = 5.0 fixed. For comparison, the figures also include a curve showing the situation in full 
LQCD with K = 0.15, NtL^ = 4^ and /3 = 5.0. The corrections up to order O (P)) in the effective 
action shift the half-filling point in the effective model from the static quark value, ^ « 1.2, a bit closer 


average Polyakov loops <L>, <L > 
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towards ^ « 1.3, the value found in full LQCD. 




Figure 19: Average sign as a function of for a two-flavor system of spatial volume V = = 4^ and with k = 0.15. 

The solid line in both figures guides the eye through the data points obtained by a Monte Carlo simulation of full two- 
flavor LQCD at inverse gauge coupling j3 = 5.0 and with Nt = 4. The other curves were obtained by a mean-field 
calculation using our improved effective action and (3.76). The left-hand figure shows results for the effective model 
for three different values of /?, with fixed Nt — 4, while in the right-hand figure, the different curves correspond to 
different values of Nt with /? = 5.0 fixed. The corrections up to order O (P)') in the effective action shift the 

half-filling point from the static quark value, /r ~ 1.2, a bit closer towards the value /r ~ 1.3, found in full LQCD. 

However, the value of the hopping parameter k = 0.15 used in Figs. 17 and 18 is already almost 
too large to be used in a calculation based on a O(jP)') truncated nearest-neighbor interaction. For 
larger values of k or lower temperatures (i.e. larger Nt), the deviation of the mean-field result from the 
Monte Carlo data becomes larger and it is no longer possible to match the data by adjusting the inverse 
gauge coupling in the effective model as was done in the right-hand part of Fig. 17. This is not sur¬ 
prising as, concerning the physics, the O{/3)) truncation of the nearest-neighbor interaction term, 
detc,d,t (l ~ implies that mesons and baryons (the low energy degrees of freedom of the theory) 

remain essentially static and cannot undergo spatial hoppings. But when lowering the temperature and/or 
the quark mass, it is precisely the hopping of these low energy degrees of freedom that becomes more im¬ 
portant. From a more technical point of view, looking at the expansion of the fermionic nearest-neighbor 
interaction term, 

detc,d,t(l - kIMc^) = 1 - Ks 

- Y (t4d.t(^C2) - 3trc,d,t(Mc2) tr c,d,t [MeJ -F 2trc,d,t(M^j) 

-F ... 

12nt 

= , (3.80) 

where the Ck{A) for a rank-n matrix A are defined by the characteristic polynomial of A, 

n 

XaW = det(lA-Fl) = 5] X'^CkiA) , (3.81) 

fc =0 

and can be obtained with the already mentioned Leverrier-Faddeev method, i.e. by starting with a matrix 
Bq{A) = 0 and the constant c„(fI) = 1, and then using the recursion 

BkiA) = ABk-i{A) + lCn-k+i{A) , Cn-k(A) = tr{ABkiA)) , k = (3.82) 

it becomes clear that the term does not necessarily give the dominant contribution to (3.80) as the 
Ck {Mc 2 ) can become rather big (especially if nt is large) and to guarantee that the K^-term dominates the 
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expansion, one would have to choose much smaller than we did. 

To incorporate spatial meson and baryon hopping into the mean-field calculation, one should include at 
least the K^-terms in the expansion of (3.80). Although this is probably feasible, already at order the 
computation by hand of the corresponding terms in the effective action becomes a delicate issue and it 
would be desirable to automate the combinatorics required to get the correct dependency on a/(/3). 

In order to demonstrate that the deviation of the mean-field result from the Monte Carlo data in Figs. 17 
and 18 is merely due to the truncation of the nearest-neighbor interactions and not caused by our mean-field 
method itself, we try to reproduce figure 6 of [6], which shows (L) and {L*} for a single flavor system with 
K = 0.01, Nt = 200 and /? = 5.7, obtained once by a Monte Carlo and once by a Complex Langevin 
simulation of the effective model derived in [6]. Our mean-field action is obtained by applying (B.43) to 
the partition function (B.40), which is of order in Ks, while for [6, Fig. 6] also terms proportional to 
have been taken into account. Flowever, as k is rather small, we should nevertheless get a result comparable 
to [6, Fig. 6] if our mean-field method works correctly. In Fig. 20, left, we show [6, Fig. 6] superimposed 
with the corresponding result from our mean-field calculation, where also for this single flavor system, we 
determined the mean-field value L within the phase-quenched system and then used reweighting to obtain 
(L) and {L*) in the phase-unquenched case. As the figure shows, our mean-field results (dotted lines) match 
remarkably well the data obtained by full Monte Carlo and complex Langevin simulations. On the right 
hand side of Fig. 20, we also compare the mean-field result obtained with the ©(k^) truncated effective 
fermion action with the corresponding result obtained with the "full" (O (P)) truncated) effective 

fermion action (B.33). As can be seen, even for the small value of k = 0.01 and large Nt = 200, the higher 
order terms in the effective fermion action lead to a clearly observable difference in the result. 
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Figure 20: The left-hand figure shows a comparison of [6, Fig. 6], showing the Polyakov loop (L) and inverse Polyakov 
loop (L*) as a function of a baryon chemical potential ps = 3p (p: quark chemical potential), obtained once by 
a Monte Carlo (MC) and once by a complex Langevin (CL) simulation of the effective model derived in [6] with 
K = 0.01, Nt = 200 and j3 = 5.7, and the corresponding result obtained with our mean-field method (dotted lines). 
As can be seen, our mean-field result coincides remarkably well with the Monte Carlo and complex Langevin data. The 
figure to the right shows again the curves that were used in the left-hand figure (dotted lines), but this time, together 
with the mean-field result obtained using the C7 (ks) truncated effective fermion action (dashed lines). As can be seen, 
even for a small value of k = 0.01 and large Nt = 200, the higher order terms in the effective fermion action lead to a 
clearly observable difference in the result. 


4 Summary 

After a short introduction to the sign problem, we have shown that the half-filling point, i.e. the value of 
the chemical potential where the fermion density assumes half its maximal value, has similar properties in 
LQCD as in e.g. the fermionic Hubbard model used in solid states physics: the sign problem almost disap¬ 
pears at half-filling and the system possesses an apparent particle-hole symmetry. We traced the origin of 
these properties by analyzing the heavy-dense fermion determinant and the corresponding strong coupling 
partition function. However, in contrast to the Hubbard model, where the lattice and the half-filling state 
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are physical, in LQCD the lattice is just a regulator that has to be removed in order to extract continuum 
physics and the state corresponding to half-filling of that lattice is therefore of no physical relevance. 

As the half-filling state in LQCD and the corresponding maximum in the average sign are just lattice arti¬ 
facts, it seems sensible that the first appearance of a minimum in the average sign (as a function of increasing 
chemical potential) indicates the point where lattice artifacts start to become dominant and measurements 
of observables taken at larger values of the chemical potential should therefore be considered with caution. 
Finally we used a mean-field calculation to obtain in an analytic way the (T, p) phase diagram for two- 
flavor LQCD (for both, isospin and quark chemical potential) in the static quark limit but at finite inverse 
gauge coupling /3, where the /3 dependency was implemented by an effective Polyakov loop action based 
on the leading terms of the character expanded gauge-field Boltzmann factor. We also showed how this 
calculation can be improved by taking into account spatial fermion hoppings. 
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A Single Site Partition Function for Mass Degenerate QCD in the 
Strong Coupling Limit 

For degenerate quark masses, i.e. Ku = Hd = k, the strong coupling single site partition function (3.18) 
with an isospin chemical potential ji = reads: 


Zs,q{fJ.,K,nt) = 


27r 277 

(^J J d0id92H{9i,62)Bet{D{0i,62; nt))'Det{D{0i,62; 

0 0 

277 277 

= (^J J deid02 H{0i,e2)\Bet{D{ei,02', fii, K,nt))\'^'^ 

0 0 

= + l){h^^+4:h^°+ 8h^+ 10h^+ I2F 

+ 22h® + 12h® + 10h'‘ + 8 P + AP + 1) 
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+ 4 (/i^ + h^) (2/)12 8^11 27/)^° + 72/)® + 115/)® + 172/)^ 

+ 193/)® + 172/)® + llhh^ + 72/)® + 27/)® + 8/) + 2) 

+ (/i® + h^) (lOh^® + 60/i“ + 20A¥° + 460/i® + 907/i® + 1232/i^ 
+ 1398/)® + 1232/)® + 907/)^ + 460/)® + 204/)® + 60/) + 10) 

+ 4(/i^ + /i®) (3/)^® + 22/)“ + 67/)^° + 172/)® + 308/)® + 450/)^ 
+ 487/1® + 450/1® + 308/i^ + 172/1® + 67/i® + 22/i + 3) 

+ 2/1® (ll/)i® + 38/)“ + 163/)^° + 386/)® + 699/)® + 974/)^ 


+ II 6 I/ 1 ® + 974/i® + 699/i^ + 386/i® + 163/i® + 38/i + ll) 


/i=(2Ke^‘)"‘ 


/i=(2Ke“f‘)"* 


(A.l) 
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The corresponding partition function with a quark chemical potential reads: 


Zsin,K,nt) 


Z7T Z7T 

J J d9id92H{0i,02) I)et{D{9i,92; /i, k, ut)) I)et{D{9i,92; /i, k, n*)) 
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ti=(2Ke^‘)"‘ 

ii=(2fce 


(A.2) 


B Higher Order Corrections to the Effective Polyakov Loop Action 

In this section we derive an effective nearest-neighbor Polyakov loop action that includes not just the 
nearest-neighbor effects coming from the gauge field action, but also those coming from the fermion deter¬ 
minant. 


We start with the single flavor partition function 


Z = 


J V[U] Det(D) , 


(B.l) 


which, using the spatial loop expansion (3.78) for the fermion determinant and the character expansion 
(3.61) for the gauge field Boltzmann factor, can be written as: 

Z = f V[U] Det(l -kT) n n n|i+E dr ariP) Xr{Up)\ ■ (B.2) 

d so {C,J p J 


We are only interested in the nearest-neighbor interactions and restrict the loop lengths sq in (B-2) to sq < 2. 
The integrand can then be written as 


Det(l — kT) 


detc.d,t(l 

<$.y) 



(B.3) 


where x, y are nearest-neighboring sites spanning a path C 2 . Equation (B.3) would capture the full spectrum 
of nearest-neighbor interactions but is unfortunately still rather complicated. We simplify the expression 
further by expanding the determinants detc,d,t (l — to order and by considering the gauge field 
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only in the fundamental representation and only along temporal plaquettes pt- With these simplifications, 
the two products in (B.3) turn into the following expression; 

{l - \ ]^|l + dfUfiP) trc(C/pj| . (B.4) 

($.y> ^ Pt ^ 

The Boltzmann factor for the desired effective action should then coincide with (B.4) up to order O Oy* (^)) 
(the single site/non-interaction terms coming from the factor Det(l — kT) in (B.2) are not considered as 
part of the effective action). 

To hnd such an action, we start by writing out tic,d,t {M^ y^ with y — x + i and * S {1, 2,3}: 


trc.d.i {Mx,y) = trc,d,4 {1 - kT)-^ {1 - kT)-^) 

nt- l ^ 

X4,y4=0 

^ ^ {Pyix4,y4) (1 + (2«;e'^)"*Pp(y4))-') 

+ (1+74) O {Pliy4,X4) (1 + (2Ke-'^)"‘pi(2/4))”')) 

• (1 +7i) 0 Uj{x,y4) 

(1 -74)0 (Ps( 2/4,X4) (1 + (2«e^)”*P,(x4))-') 

+ (1 + 74) «) (4(^4,y4) (1 + (2Ke-'^)"‘pi(x4))”')) | . 

(B.5) 


After carrying out the trace in Dirac space this yields: 

nt- l / 

trc,d7(M$.y) =2 ^ 

X4,y4=0 ' 

■ tic{Ui{x,X4) Py{x4,y4) (l + (2«;e'")"*Pp(j/4)) ^ 

• Uj{x,y4)Ps{y4,X4) (1 + (2Ke'")"‘P5;(a:4))”^) 

— (2 (2 ^_2^^e(a;4>y4)+e(y4<a;4) 


• trc(Pi(a;,a:4) Pp(a:4,y4) (l + (2 H;e^)”*Pp(t/4)) 

• U^^{x,y4) Pl{x4,y4) (1 + (2Ke"^)”‘P^(x4))”^) 


• trc(Pi(a;,X 4 ) Pg(?/ 4 ,a; 4 ) (l + (2/ce '")"'‘pi(?/ 4 )) ^ 

• Uj{x,y4) Pxiy4,X4) (l + (2Ke^)”‘Ps;(a:4)) ^) 

+ (2 K ^2 K ^_^'^e{y4>X4)+e{xi>yi) 

■ tTc{Ui{x,X4) pl{y4,X4) (1 + (2/ce“'")"'‘Pj(j/4)) ^ 

•P;(x,t/4)Pi(x4,2/4)(l + (2Ke-'^)”‘P-^(x4))”')) . (B.6) 
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Integrating over the spatial links in (B. 6 ) at order leads to the constraint j /4 = X 4 (see Fig. 21) 

as for SU(3) integrals, we have; 


JdUUa,b=0 , J dU Ua,b Uld = i 6a,d Sb,c • (B. 7 ) 

According to the dehnition of Px{x4, yy) given below eq. ( 3 . 29 ), we then have P^{x4, X4) = 1 and hnd 
therefore: 

(tre((l + ( 2 Ke^)"‘Pj^)-')-trc((l + ( 2 Ace-^)”*Pj)-^)) + 0 (a/(/ 3 )) , (B. 8 ) 
which is the same as eq. ( 2 . 28 ) of [ 6 ], as can be seen by noting that; 

tr,(( 2 «e'^)"‘P 5 (l + ( 2 «e^)”*Ps)-') = 3 -tr,((l + ( 2 A«e^)"‘Ps)-') , (B. 9 ) 

and 

tre(( 2 Ace-^)”‘P-^(l + ( 2 Ke-'^)”‘P-^)”^) = 3 - tre((l + ( 2 Ace-'^)”*Pj)”^) . (B.IO) 






Figure 21: The figure shows the four types of diagrams corresponding to the four terms in (B.6) that give non-vanishing 
contrihutions after integrating out the spatial links in the fundamental representation (left-hand sides of the little black 
arrows), together with the corresponding configuration after the spatial links have been integrated out (righ-hand sides 
of the little black arrows). 


To obtain non-zero contributions from terms in (B. 6 ) for which j /4 ^ X 4 , we have to take plaquette 
terms into account. To see in more detail how this works, let us as an example consider the top left diagram 
in Fig. 22. It shows a contributing diagram to the first term inside the large bracket in the sum over X 4 , j /4 
in eq. (B. 6 ) for the case where y 4 — X 4 = 2. If we would integrate out the spatial links in the corresponding 
term in (B. 6 ) without taking into account plaquette terms, the result would just vanish due to the left-hand 
identity in (B.7). But allowing for plaquette terms coming from the non-trivial gauge-field part in (B.4), we 
can form the two left-most diagrams on the second and third row of Fig. 22, and get for example for the 
one on the second row, a term like 


- (3a/(^) tre(P^J) (3a/(^) trc(P^J) 

(2 %ic{Ui{x,X4) Py{x4,X4 -F 2) (1 -F {2 Py{x4 + 2)) ^ 

U^{x,X4 + 2) P^{x4 + 2,X4) {1 + {2Ke>")'"*P^{x4)) , 


where 

and 


= U4{x,X4)Ui{x,X4 + l)Ul{y,X4)Ul{x,X4) 

Pp 2 = U4(x,X4 + 1) Ui{x,X4 + 2) p|(?/,a ;4 + 1) U]{x,X4 + 1) . 


(B.ll) 

(B.12) 

(B.13) 
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Integrating over the three spatial links Ui{x, X4), Ui{x, X4 + 1) and Ui{x, X4 + 2) then yields 


-2(3ay{^)fl(2rceT‘ 

^^c{yl{y,X4 + 1) Ul{y,X4)Py{x4,X4 + 2) (1 + (2Ke'')"*Py(a:4 + 2))"^^ 

iT:c{UA{x,X4)U4:{x,X4: + l)Px{xi + 2,X4){l + {2Ke^"Y‘^Px{xi)) , (B.14) 

where, as per definition (see Sec. 3.2.3), 

Py{x4,X4 + 2) = U4{y,X4)U4{y,X4 + l) (B.15) 

and 

Ps{x 4 + 2, X 4 ) = U 4 {y, X 4 + 2 ) ... 174(1/, rit - 1) U 4 iy, 0) 1/4(2/, 2^4 - 1) , (B.16) 

we can simplify 

17|(2/,a:4 + 1) 174(2/,a;4) 7’y(a;4,2;4 + 2) = 1, (B.17) 

and 

1/4(2;, 2:4) 1/4(5:, X4 + 1) Pj(a:4 + 2, a;4) = ^5(2:4) , (B.18) 

to get either the identity or a closed Polyakov loop. With this, we finally obtain for eq. (B.14): 

- trc((l + (2^te'^)"‘Py)”^) tre((2Ke^)”‘Ps(l + (2Ke^)”‘Ps)”i) ■ (B.19) 

In a similar way we can proceed with the other terms in (B.6) up to order 0(^k^ 0/*~^(/3)) and find, 
after grouping them according to the power of their factors of af{/3): 

-I E {(a/(/3))“°"^"^-"^’”‘^(trc((l + (2«e'^)"*P5)-') tr,((2r^e^)"‘P, (l + (2r^e^)”*P,)-') 

+ tr,((2Ke-'^)"‘Pj (1 + (2Ke-^)”*Pj)-^) tr,((l + (2Ke-^)”*Pj)-') 

+ tr,((l + (2Ke^)"‘Pj^)”^) tr,((l + (2rce-^)”‘P-^)”^) 

+ tre((2Ke-'^)”‘ Pj (l + (2 «;e-'^)”‘Pj)”^) 

trc((2Ke'")"'* Pj (1 + (2Ke'")"*P5=) 

+ (^tr,((2Ke'^)"* Py (l + {2Ke^rPy)~^) tr,((l + (2rce^)”‘P5)”^) 

+ tr,((l + (2r^e-^)"*Pj)-') tr,((2r^e-^)”* P^ (l + (2r^e-^)”*Pj)-') 

+ tr,((l + (2K:e-'^)”‘pi)”') tr,((l + (2 K:e'')"‘P5)”') 

+ (2rc)-"“°‘^(""-*'"’”‘hr,((2Ke'^)”‘ Py (l + (2 Ke^)”*P,,)”^) 

tr,((2Ke-'^)"‘ P-^ (1 + (2/ce-'^)”‘P-^)”^)^ | , (B.20) 

Note that all the traces in (B.20) are time-independent. The time-dependency of the individual terms in the 
sum is only in the different exponents of af{j3), (2 and (2 k)~‘^ which form essentially simple geometric 
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sequences. The (finite) sum can therefore easily be carried out and yields: 


2 f nt(af(/3) - af{l3)) 
'3 1 l-a/(/3) 




tre((l + (2Ace^)”‘P5) ") 

trc((2Ke^)”‘ (1 + (2Ke^)”*P5;) 

+ tr,((2Ke^)"‘ Py (1 + (2Ke'")"‘Py)”^) 

tre((l + (2Ke^)"‘Ps)”^) 

+ tr,((l + (2Ke-'^)”‘Pj)-') 

tr,{{2ne-^r4{l + {2ne-^r4)-^) 

+ tr,((2Ke-^)”*Pj (1 + (2Ke-'^)”*Pj)”^) 

tre((l + (2Ace-'^)”*Pj)”^)^ 


nt{{2K) af{j3)-{2K) "'*af{l3)) 
1 - (2K)^a/(/3) 




nt{{2Kf^^afi^) - i2Kfaf{P)) 
(2k)^ - ttfiP) 


trc((l + (2Ke^)”‘Pjj 

tr,((l + (2«e-^)”*P,^)-') 

+ trc((l + (2Ke“'")"*Pj) 

trc((l + (2Ke^)”‘Pj)”^)^ 

tr,(P5(l + (2Ke'^)”‘P5)-') 

tr,(pt(l + (2«e-'^)”‘P,^)-') 
+ tr,(pJ(l + (2K:e-'^)"‘Pj)-') 

tre(Ps(l + (2Ke'^)"‘P5) 


(B.21) 


For X 4 ^ j/ 4 , there is a third class of diagrams, which are always of order a^* (/3), depicted in the last 
row of Fig. 22. If we take as an example again the left-most diagram on that row, we obtain by integrating 
out the spatial links, using the identity 


1 

d U Uai^bi Ua 2 ,b 2 ^ ^Oi 03 C&i 63 63 i 


(B.22) 


a term of the following form: 


- 2(3a/(/3))”‘ (2 k e^) 




nt -2 


e/i I2 h J2 J3 M2 M3 iNi N2 N3 


Py{xA, 2 / 4 ) j, Ml Pi (a^4, y4)N, h ( 2 ^ 4 , X4)j^ P^[y4, Xi)^^ 

{Py{x4,y4){l + (2Ke'')”‘Py(2/4))”\^jv^3 {Ps:{y4,X4){l + (2Ke"'')"‘Ps(a;4))’ 

By going to maximal temporal gauge, such that 


’N 3 I 3 


Px{x4,y4) = 


if X4 < 2/4 


Px ) else 


. (B.23) 


(B.24) 
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we can simplify: 


J 2 Jz ^M\ M 2 M 3 Pyip^A-) 2 / 4 ) Py ( 2 / 4 ? ^ 4 ) J 2 M 2 

{ f^c{^^) *^J 3 M 3 Px^M^ J 3 : if X 4 > y 4 

trc(pj) ^J 3 M 3 - M 3 J 3 . if 2^4 < y4 , (B.25) 

2 <5j3 M 3 , if a;4 = j /4 

and (B.23) becomes for X 4 > 1 / 4 : 

rLf / n\ 

- ( 2 «e'^)"*(tr,(P 5 ) tr,(P 5 (l + (2«e^)"‘Pj^)-') - tr,(P|(l + (2«e^rPj^)-')) 

(tre(pf) tr,((l + (2Ke'^)”*Ps)-') -frc(4(l + (2'«e'^)”‘P5)”')) , (B.26) 

and for X 4 < y 4 : 

(^TLt / 0\ 

- ( 2 «e'^)"‘(tr 4 pJ) tr,((l + Py)-^) - tr 4 pj(l + {2Ke^rPy)-^)) 

(tr,(P 5 ) tr,(Ps(l + ( 2 Ate'^)"‘P 5 )"') - tr,(P|(l + ( 2 Ate'^r‘P 5 =)-')) . (B.27) 

However, by using that for P G SU(3), we have 

pt = p 2 _ tr4P)P + 1 trpf , (B.28) 

we see that (B.26) and (B.27) are actually the same, and can be written as, 

- ( 2 .«e'^)"‘(tr,(P 5 ) tr,(P 5 (l + (2.^e^)"‘P 5 )-') - tr,(P|(l + (2Ken^‘Py)-)) 

(tr,(P 5 ) tr,(P:j(l + (2Ke'^r‘P:,)~^) - tr,(P|(l + ( 2 Ate'^)"‘P 5 )-')) . (B.29) 

The terms corresponding to the remaining diagrams in the last row of Fig. 22 can be computed in a com¬ 
pletely analogous way. Collecting all the resulting terms and summing over X 4 ^ j /4 then yields; 

af‘(/3)nt I 

(nt - l)(2«e'^)"‘ (tr,(P 5 ) tr,(Pg(l + (2Ke'^r^Py)~^) - tr,(P^^(l + (2 Ken"^Py)~^)^ 

■ (tr^P^) tr,(P 5 (l + ( 2 Ke'^pP^)-) -tr 4 P^^(l + (2 Ke'^pP^)-)) 

+ -(i^c(Py) tr 4 Py(l + ( 2 Ke'^rPy)-')-tr,(Py^(l + {2 Py)-")) 

1 — I ‘2iHj 

■ (trc(Pi) tr,(pf(l + (2K;e-'^)"‘pf)-') - tr,(Pf (l + (2«e-^)"‘pf)-')) 

• (tr,(Ps) tr,(P 4 l + ( 2 «e^)"*Ps)-') - trc(P|(l + {2 n P^)-")) 

+ {nt - l)( 2 «;e-'^)"‘ (tr^Pj) trc(Pj (l + (2 Ke-^)”*Pj)”^) - trc(pi^(l + (2 Ke-'^)”‘Pj)”^)) 
■(trc(Pi)tr,(pi(l + (2^e-'^)”‘pi)-Vt^c(Pf(l + (2«e-^)"‘pi)-'))|. (B.30) 
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Putting everything together, the effective nearest-neighbor Polyakov loop action, whose Boltzmann 
factor coincides with the full nearest-neighbor interaction term in (B.4) up to order a^*{/3)), can be 
written as 

- Seff{nt,K,l3,fj,) = -Sg^eff{rit,/3) - Sf^eff{nt,K,P,fi), 

where 

- ^ (trc{Ps) trc(P’J) -f trc(P^) trc(P5)) 

{x,v) 

is the gauge part, already known from (3.62), and 

- Sf,eff{nt,K,f3,fl) = 

2kI nt 


(B.31) 

(B.32) 


^ |('tr,((l + ( 2 /ce'^)”‘P 5 ) ^)-tr,((l + ( 2 Ke-^)"‘pi) 

(s.y) ^ ^ 

• ^trc((l -f {2Ke^)'^*Py)~^) - trc((l -f (2 

(trc((l + (2Ace^)"‘P^)-')(3-tr,((l + (2 Ke'^)”*P,)-')) 
+ (3 - tr,((l + (2Ke^)”‘P5)”^)) tr,((l + (2Ke'^)”‘P 5 )”^) 


+ tr,((l + (2«e-'^)"‘Pj) )(3-tr,((l + (2«e-^)"‘pT) )) 

+ (3- tr,((l + (2«:e-'^)”‘Pj)-')) tr,((l + (2/re-'^)”‘P,^)-')) 

(2^) ^^V tr,((l + (2«e^)"‘P,)-') tr,((l + (2«e-^)”*Pj)-^) 

l-(2K)a/(/3) V 

+ trc((l -f ( 2 «;e"'")"*Pj)”^) trc((l -f (2 Ke'")”*P 5 ;)”^)^ 


_ i2n) a^iP) (2^) «/(^) ^r,(P,(l + (2 P,)-^) tr,(pj(l + (2«e-'^)"‘Pj)-') 

( 2 k) -afiP) \ 

+ tr,(pJ(l + ( 2 Ke-'^)”‘pi)-')tr,(P^(l + ( 2 Ke^)”‘P,)-')^ 

(tr^P^) tr,(P^(l + (2Ke^)"‘P^)-') - tr,(P|(l + (2Ke^)"‘P^)-')) 

• (tr,(P^)tr,(Ps(l + ( 2 Ke^)”‘Ps)”Vtrc(J^I(l + ( 2 Ke'^)”‘Ps)’ 


3a/(/3) ((2K) (2k) ) ^ (2Ke^)"‘P,)-^) - tr,(P|(l + (2Ke^)"‘Pg)-')) 

4 (l — [2k) ) ^ ^ 


(tr,(pt) tr,(pt(l + ( 2 Ke-'^)"‘P,^)-') - tr,(pf (l + ( 2 Ke-^)"‘P,^)-')) 


^ (t^c(P^)tr,(pt(l + ( 2 Ke-")"P^)-Vtrc(pj^(l + ( 2 Ke-")"pt)-^)) 

• (tr,(Ps)tr,(P 5 (l + ( 2 Ke^)"*Ps)-') -trc(P|(l + ( 2 Ke^)”*Ps)-')) 
3(n,-l)a/(/?)(2Ke ^ (2Ke-^)”‘Pj)-') - tr,(Pf (l + (2Ke-'^)”*Pj)-')) 


(tr,(pt) tr,(pt(l + ( 2 Ke-'^)”‘P, 0 -V trc(Pr(l + ( 2 «e-^)"‘Pi)-')) | 


(B.33) 


is the fermionic part. Note that the k^Oj* (/3) contribution and the k'^ contributions are already contained in 
the general effective action derived in [6], [11] and [12]. 
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Figure 22: On top of the figure, we show four diagrams corresponding to the four different terms in (B. 6 ) for a particular 
choice of 2:4 and 2/4 (with y^ — Xi = 2). Below each diagram, the three possibilities of attaching fundamental plaquettes 
(coming from the gauge action) to the diagrams, such that they survive the process of integrating out the spatial links, 
are shown together with the resulting loop configuration after that integration. Solid lines correspond to terms like 
(B.36) and (B.37), dashed lines to terms like (B.34) and (B.35) and the mixed solid and double solid lines to terms like 
(B.38) and (B.39). 


By using 


and 


trc((l + (2Ke = 


3 + 2(2Ke^)”‘ tic{P) + tr^P^') 

1 + (2k trc(P) + (2 k trc(P'l') + (2Ke^)^”‘ 

3 + 2(2Ke-'^)”‘ trc(Pl') + (2Ke-'^)^”‘ trc(P) 


(B.34) 


1 + (2Ke“^)"* trc(P'l') + (2Ke“'^)^"‘ ^^c{P) + (2Ke“^)^”* 

(B.35) 
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as well as 


(B.36) 


trc(P(l + ( 2 Ke^)”‘P) 


-1\ 


3(2Ke^)^”‘ +trc(P) +2(2Ke'")”‘ 

1 + ( 2 k e^)"* trc(P) + ( 2 k trc(P'l') + ( 2 k e'^) 


and 


tr (pt(l + (2Ke--r‘Ptr^) = 3(2Ke-^)^"-+tr.(pt)+2(2Ke-^rtr4P) 

^ ^ ’ l + (2Ke-/^)”*tr,(Pt)+(2Ke-'^f”‘tr,(P)+ (2Ke-A‘)'"‘ ’ 

(B.37) 

and also 


(tre(P) tr,(P(l + (2Ke^rP)”^) - tr,(p2(l + (2Ke^)”‘P)”^)) 

_ 2 trc(Pl') + (2Ke'^)"‘(3 + 2 ( 2 Ke^)”‘tre(P) +trc(P)trc(pt)) 

1 + (2Ke^^)”‘ tr<,(P) + (2Ke^)^”* tr^Pt) + (2Ke^)^”* 


(B.38) 


and 


(B.39) 


(trc(P^) trc(P^(l + ( 2 Ke"'")”‘P^) ^) - trc(P^^(l + ( 2 Ke"'")”‘P^) ^)) 

_ 2 trc(P) + (2Ke-'^)”*(3 + 2(2Ke-^)”*trc(p1') + tr^pt) tr^P)) 

1 + ( 2 Ke“'^)”‘ trc(P'l') + ( 2 Ke“'')^”‘ trc(P) + ( 2 Ke“'^)^"* 

equation (B.33) can be written completely in terms of Polyakov loops Li = trc(Pj.) and complex conju¬ 
gate (or inverse) Polyakov loops L* = trc(Pj ) at neighboring spatial sites Xi, Xj, i ^ j. 

The partition function for the effective single flavor theory can now be written as 


Z,ff{nt,K, P,fi) = J V[P] I (jl det2(l + (2Ke'^)"‘P5) det^ (l + (2 k e-^)”*Pj)) 

where the correspondence between the terms in the integrand of (B.40) and the terms in the exact partition 
function (B.l) is given by: 

Det(P) - (j|det^(l + (2Ke'")”*Ps)det^(l + (2Ke-'')"‘pi)j (B 41 ) 

^ X ' 

and of course: 

e-'S's ^ Q-Sg.,ffint,P) _ (B.42) 

The generalization to Nf flavors is, at the present expansion order, obtained by simply adding appropriate 
factors of (B.41) to the integrand of (B.40). 


(B.40) 


The single site mean-field action corresponding to (B.31), which was used to generate the figures 17, 18 
and 19 in Sec. 3.6, is obtained as usual, with the small complication that, as the action is no longer bi-linear 
in the Polyakov loop variables but a rational function, the usual approximation procedure: 

1. write L, = L + SLi, L*=L*+ SL* in Seff, 

2. keep terms only up to linear order in the perturbations 5Li, SL* 

3. substitute back SLi = Li — L, SL* = L* — L*, 

4. set Li = La, L* = La for all i (the subscript a means "active") and divide by the system size. 
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should be interpreted as coming from a leading order Taylor expansion of S^f / around the state where the 
Polyakov loops and inverse Polyakov loops at each spatial site assume their mean-field values, i.e.; 


E{(La-L) 

i ^ 


, Ly Ly) 


SmfiL^.LlL) = 


Li^L 



where V = n^nyUz is the spatial volume and the leading constant piece has been dropped. Note also 
that the static-quark part of the fermion determinant, i.e. the product of single site static-quark fermion 
determinants in (B.41) is not part of the effective action and therefore not expanded during the mean-field 
approximation. 
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